Quark mass variation constraints from Big Bang nucleosynthesis 
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We study the impact on the primordial abundances of light elements created by a variation of 
the quark masses at the time of Big Bang nucleosynthesis (BBN). In order to navigate through the 
particle and nuclear physics required to connect quark masses to binding energies and reaction rates 
in a model-independent way, we use lattice QCD data and a hierarchy of effective field theories. 
We find that the measured 4 He abundances put a bound of —1% < 8m q /m q < 0.7% on a possible 
variation of quark masses. The effect of quark mass variations on the deuterium abundances can be 
largely compensated by changes of the baryon-to-photon ratio n. Including bounds on the variation 
of r\ coming from WMAP results and adding some additional assumptions further narrows the range 
of allowed values of 5m q /m q . 



I. INTRODUCTION 

In theories of physics beyond the standard model the 
standard model parameters appear not as fundamental 
constants but as derived quantities. In many of those 
theories the possibility then arises that the values of the 
standard model "constants" can vary over time [I . It is 
then important to understand which constraints the suc- 
cesses of standard cosmology - which assumes time inde- 
pendent constants - poses on this purported time varia- 
tion. A natural place to look for a strong sensitivity to 
a variation of fundamental constants is Big Bang nucle- 
osynthesis (BBN) since it satisfies two important criteria. 
First, BBN happened very early in the universe's history, 
mostly when the universe was between 3 seconds and 3 
minutes old. Second, not only is standard BBN under- 
stood at a few percent level but it is also very sensitive to 
microscopic parameters such as nuclear binding energies 
and reaction rates that are, themselves, very sensitive to 
certain standard model constants. It is no surprise then 
that BBN has been used in the past to study the vari- 
ation of fundamental constants [2]- The purpose of the 
present paper is to explore the BBN constraints on the 
variation of the masses of the two lightest quarks, m u 
and mj. 

The binding of nucleons into light nuclei during BBN 
proceeds through a number of reactions, some of which 
are in equilibrium with the expansion of the universe and 
some that are not. After weak reactions like p + e~ -H> 
n+v are no longer in equilibrium (i.e. weak freezeout) , the 
ratio of neutron to protons decreases due to neutron j3- 
decay. If the formation of light nuclei occurred in equilib- 
rium, the most bound nuclei (among the light ones this is 
4 He ) would form earlier and more abundantly. The for- 
mation of 4 He can, however, only occur after 2 H , 3 He and 
3 H have been formed, since multinucleon fusion reactions 
are essentially impossible at the relatively low densities 
prevalent during BBN. Their number is small on the ac- 
count that their binding energies are small and it is not 



energetically favorable for them to form until the temper- 
ature is low enough to be comparable to their binding en- 
ergies. Thus, the beginning of nucleosynthesis is delayed 
by the shallowness of the deuteron binding energy, the so- 
called deuterium bottleneck. Since this shallowness is a 
product of delicate cancellations between kinetic and po- 
tential energies, the binding of the deuteron is an obvious 
place where a small change in quark masses can signifi- 
cantly alter the primordial abundances. Notice that the 
rate for the reaction n + p O d + 7 is not small; it is suf- 
ficient to keep the deuteron number in thermal/chemical 
equilibrium. It is the equilibrium deuteron number that 
is too small for them to collide and be assembled in larger 
nuclei. After the deuteron number grows enough, the re- 
actions leading to the formation of 4 He proceed quickly 
and essentially all the neutrons present in the beginning 
of BBN are assembled into 4 He nuclei. The timing where 
this assembly starts (determined, among other things, by 
the deuteron binding) is crucial as the neutron numbers 
are decreasing due to neutron /3-decay. Small amounts 
of 2 H , 3 He and 3 H are left out of this process. Their 
numbers depend critically on chemical non-equilibrium 
physics and the rates of the reactions, including the ini- 
tial n + p — > d + 7 reaction. Current observation is not 
useful in measuring reliably the primordial abundance of 
3 He and 3 H . However, the abundance of 2 H , and es- 
pecially 4 He , are well measured and put a significant 
constraint on any change of the standard BBN scenario. 

A number of authors have previously considered the 
effect of quark mass variations on the BBN predicted 
abundances [3TU5| . The main difficulty to be surmounted 
is that the quark mass dependence of binding energies, 
reaction cross sections and decay rates that are input 
to BBN models are difficult to determine. For instance, 
modern nuclear potentials can describe very well nucleon- 
nucleon phase shifts. They can also be used to compute 
binding energies with enough precision (with the help of 
phenomenologically motivated three-nucleon forces fit to 
some observables) and cross section for few-nucleon reac- 
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tions. These potentials are, however, tuned to data ob- 
tained from experiment where the quark mass has its cur- 
rent value. What is usually done in estimating the effect 
of quark mass variation is to change the parameters in 
these models where this dependence is easy to track. For 
instance, the range of nuclear forces, given by l/m„-can 
be changed through the relation m\ ~ m q . But the long 
distance part of the potential, sensitive to this range, is 
actually a small part of the nucleon-nucleon interaction. 
The medium and short range parts also have a quark 
mass dependence and, while this dependence is likely to 
be milder, its effect on the overall nucleon-nucleon inter- 
actions is still large due to fine-tuned cancellations that 
are responsible for, among other things, the shallowness 
of the deuteron. In this paper we avoid, as much as 
possible, model- dependent approximations of the proper- 
ties of nucleons and its nuclear forces, relying solely on 
the symmetries of QCD and its connection to nuclear 
physics through more general arguments. In particular 
we use effective field theories (supplemented with lattice 
QCD data) to connect the change in quark masses to the 
inputs used in BBN simulations. 



A. Effective Field Theories 

At momentum scales Q below Aqcd ~ 1 GeV, the 
relevant degrees of freedom in QCD are hadrons, not 
quarks and gluons. Effective field theories (EFTs) for 
this momentum range (i.e. chiral perturbation theory) 
were developed for the meson, one and many-nucleon 
sectors. They are able to predict physical observables 
as an expansion in the small parameters Q/ Aqcd an d 
TOjr / Aqcd, Taking as inputs a few "low energy con- 
stants" (LECs), like pion decay constants and the nu- 
cleon mass in the chiral (m^ = 0) limit. These LECs, 
in turn, are determined from analyses of experimental 
results. Effective field theories predict, for example, the 
dependence of nucleon masses on the value of the quark 
masses. This particular change, however, is very small 
and can be neglected, except for its effect on the phase 
space for the neutron decay and related weak process (see 
below). Lattice QCD calculations reinforce the belief in a 
small quark mass dependence of nucleon masses [THl HZ] • 
On the other hand, chiral perturbation theory for few- 
nucleon systems (referred to here as %EFT) is in a less 
developed phase. First, there are conceptual issues that 
preclude a reliable prediction of the quark mass depen- 
dence of few-nucleon observables [IS]- Second, it has not 
been used extensively in multi-nucleon systems and their 
reactions involving photons. To bypass this difficulty we 
use a low-energy effective theory where all particles, such 
as the pions, have been integrated out with, leaving only 
the nucleon, photon, and neutrino degrees of freedom. 
Known as "pionless EFT", the momentum scales Q rel- 
evant to this effective theory are much smaller than the 
pion mass m w . This theory can make non-trivial pre- 



dictions because the states pertinent to BBN ( 2 H , 3 He , 
4 He ) are loosely bounded and the typical momenta Q 
of their constituents are significantly below ~ and 
therefore within the regime of validity of the pionless 
EFT [10]. The pionless EFT is very successful in pre- 
dicting observables in the three-nucleon sector and there 
is indication that the same is true in the four-nucleon 
sector [HO HO]. Since the a-particle is the most bound 
of s-wave nuclei, its successful description in the pion- 
less EFT might indicate that the theory can be useful in 
studying larger nuclei. Since the pionless EFT makes no 
use of the QCD chiral symmetry, it cannot directly pre- 
dict the quark mass dependence of observables. The pa- 
rameters of the pionless EFT, at the lowest orders in the 
low energy expansion, are the threshold nucleon-nucleon 
scattering parameters (e.g. scattering lengths, effective 
ranges, etc. . . ). These few parameters have been stud- 
ied using %EFT and we can use them to predict their 
variation with quark masses. In addition, some lattice 
QCD results confirm and reinforce the xEFT predictions 
for scattering length dependences on quark masses. We 
use these %EFT results as input parameters for the pi- 
onless EFT. This allows us to obtain estimates for the 
quark mass dependence of nuclear properties relevant to 
BBN. We will then use this information in combination 
with a standard BBN code to compute the light elements 
abundances in order to constrain the values of the quark 
masses during the Universe's first minutes. Our strategy 
of combining these two types of effective theories is sum- 
marized in Fig.[T] We will now describe the stages of our 
calculation. 



B. Scattering length dependence on quark masses 

Different versions of %EFT have been used by different 
authors to study the quark mass variation of the nucleon- 
nucleon S*-wave scattering lengths. The results depend 
on the spin channel. In the spin singlet 1 5ochannel and 
at leading order (LO) on the m^/ Aqcd expansion, the 
calculation of the quark mass dependence of the scatter- 
ing length in the version of xEFT used in [21] requires 
as inputs the chiral limit values of the axial charge of 
the nucleon gA, the decay constant of the pion /„., the 
nucleon mass M, the pion mass m v and the coefficient 
of a two-nucleon contact term C° fitted to the physical 
scattering length. Only the value of these quantities at 
the physical value of quark masses is precisely known, 
but the difference between them and their chiral limit 
values is a higher order effect that can be neglected in 
a next-to-leading order (NLO) calculation. At NLO a 
new constant D 2 S appears (which is the coefficient of a 
two-nucleon operator with no derivatives but one quark 
mass insertion) as well as other constants contributing 
to the quark mass dependence of f v ,gA an d M. The 
value of Dg is difficult to disentangle from C° as both 
contribute equally to nucleon-nucleon scattering at the 
physical value of the quark masses. They give, however, 



3 



binding energies 

§Ai /71- j M , ■ • ■ ==> a2 , at , • • • . . . _^ abundances 

XBFT pionless EFT reactivities BBN code 

FIG. 1: The strategy used to determine quark mass dependence of BBN abundances. At the far left, empirically determined 
LECs are used to constrain xEFT, allowing predictions of nuclear observables and determination of their quark mass dependence. 
This theory is in turn matched onto the pionless EFT, where subsequent calculations of binding energies and reactivities relevant 
to BBN are used as input into BBN codes. 



different extrapolations to other values of quark masses. 
They can be disentangled only through a study of pro- 
cesses like deuteron-pion scattering or by the use of lat- 
tice QCD data (see below). The strategy used in deal- 
ing with the lack of knowledge of the value of D 2 is to 
estimate it using naive dimensional analysis arguments. 



In Ref. 



Di was constrained by requiring its abso- 



lute value not to be too much larger than \Cg\ while 
in Ref. [22j D| itself was required to be of natural size. 
Fortunately, the difference in the power counting schemes 
used in [5T] and [55] has little impact on the dependence 
of the scattering length on the quark masses and the dis- 
crepancy between them can be explained by the different 
assumptions about the reasonable range of values for D 2 . 
We will use the calculation described in [55] as those au- 
thors computed the quark-mass dependence of both the 
deuteron binding energy and nucleon-nucleon scattering 
lengths, since the deuteron binding energy is one of the 
most important ingredients in the BBN calculation. 

For a small variation of the quark mass we can read 
off figure If in Ref. [22 the slope (we use the more con- 
servative estimate where the change of the axial constant 
qa with quark masses, parametrized by d\Q is included): 



m q 8B2 SB2 
B2 5m„ 2B2 Sm n 



-f (-0.085 ±0.027), (1.1) 
its 2 



where m q is the average mass of the up and down quarks 



and we made use of the relation m 2 



Similarly, we 



use figure 12 in [55] to extract the variation of the spin 
singlet 1 channel scattering length to find 



Sa s 2m q Sa s 



(-1.4 ± 1.4) 



fm 
MeV 



(1.2) 



Notice that a vanishing a s variation is consistent with 
these extrapolations, a feature also seen in the extrapo- 
lation in |21j . If a s were the only parameter determining 
the change of abundances due to varying quark masses, 
BBN would impose no constraint on possible quark mass 
variations. 

Fully dynamical lattice QCD calculations of the 
nucleon-nucleon scattering lengths have appeared in the 
last few years. They are still performed at higher values 
of quark masses, too high for the effective theory ap- 
proach to be valid, so they are of limited value for our 
purposes. Despite that, an attempt was made in [23] to 
use xEFT to find the quark mass dependence of scatter- 



ing lengths by interpolating the lowest pion mass lattice 
data and the known experimental value of the scattering 
lengths at the physical point. At this point in time, the 
deuteron binding energy has not been measured from lat- 
tice QCD. However, it is related, at leading order in the 
effective theory, to the triplet scattering length that is 
measured. Using the extrapolation in 1231 and the lead- 
ing order relation B 2 = I /(Ma 2 



if) we find 



TTlq dB 2 

B2 drriir 



= -0.14 ±0.13, 



(1.3) 



in agreement with Eq. ( |1.1[ ). In the extrapolation done 
in |23] another branch of allowed values of dB2 j dm^ ap- 
pears. This additional band is excluded from the purely 
EFT extrapolations in [5T] and [55] and will be disre- 
garded in this paper. 

The allowed values for the a s quark mass dependence 
extracted from the extrapolation in [23 , namely 



(-0.75 ±1.0) fH) 



drrin 



MeV 



(1.4) 



are consistent with the ones above but are too loose to 
add any relevant constraint. 

The remaining inputs of the pionless EFT, like three- 
nucleon interaction parameters, effective ranges, nucleon 
magnetic moments, etc. . . , are not fine-tuned and there- 
fore vary much less drastically with the quark masses. 
Their contribution to the overall fusion cross sections is 
also suppressed compared to B2 and a s . In the present 
paper we will take them to be independent of the quark 
masses. 



C. Binding energies, reactivities and lifetimes 

We have used the pionless EFT to estimate the quark 
mass variation of four quantities: the binding energies of 
the deuteron, 3 H , 3 He , 4 He and the reactivity of the pro- 
cess n + p — > d + 7. Similar calculations were carried out 
for 3 H in the context of infrared limit cycles in Refs. [541 - 
126] . The binding energies of larger nuclei, like 7 Li , are 
important only for the abundances for these larger nuclei. 
As it is not presently possible to have a reliable estimate 
on the quark mass variation of these binding energies, 
we keep them fixed and concentrate on the abundances 
for the lighter nuclei 2 H and 4 He , confident that they 
will not be significantly affected by the binding of A > 4 
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nuclei. We also only include the variation of the reac- 
tivity of proton-neutron capture as this is the reaction 
that initiates BBN and is more likely to have an impact 
on abundances (but, as we will see below, this impact is 
minimal). The binding energy of the deuteron is given 
by Eq. (Ll). 

The calculation of three-nucleon and four-nucleon 
properties in the pionless EFT requires as inputs the sin- 
glet and triplet scattering lengths as well as one three- 
body observable, usually taken to be the triton bind- 
ing energy. This binding energy can be traded by the 
value of a three-body force counterterm. The three-body 
force is also not fine tuned and will therefore show only a 
weak dependence on the quark masses that we will con- 
sequently neglect. Changing the two-body input while 
keeping the three-body counterterm fixed provides then 
the scattering length dependence of the three-nucleon 
system. In other words, the binding energies of the 3 He , 
3 H and 4 He nuclei are estimated by 
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explicit formula from [27] that the scaling with B 2 is 
indeed very well satisfied but that the scaling with a 2 
does not work as well. 

Finally, we discuss how quark mass changes affect the 
neutron lifetime as well as the rates of other one-baryon 
weak reactions such as p + e~ <-> n + v. This effect 
arises from a modified value of the axial charge and 
the neutron and proton masses, which in turn dictate the 
allowed kinematic phase space for these weak reactions. 
In fact, the neutron width is given by |28) 



(G F cos 9 c f 
2^3 



ml(l + 3g 2 A )f 



A 



(1.8) 



where A = m n — m p and m e are the mass splitting be- 
tween neutron and proton and the electron mass, respec- 
tively, g A ~ 1.26 is the nucleon axial decay constant, 
Gf the Fermi constant and 8 C the Cabibbo angle. The 
function /(A/m e ) is 



m q dBi 



Bi dm q Bi \dm q da 



q ( da s dBi dB2 dBi 
dm q dB 2 



(1.5) f(w ) 



where Bi stands for the binding energy of one of 3 He , 
3 H or 4 He . The values of the derivatives appearing in 



Eq. (1.5) were computed using the pionless EFT: 



dBz 



B 3 da s 
a s dBi 
B4 da s 



= 0.12, 



= 0.037, 



B 2 dB A 
B~ 3 dB~ 3 
B 2 dB A 
IhdBz 



= 1.41, 



= 0.74, 



(1.6) 
(1.7) 



where B4 is the 4 Hc binding and B 3 the 3 H or 
3 He binding energy. The weak dependence on a s is easily 
understood when one notices that the typical momenta 
in these bound states is of order \fMB~i, which is much 
larger than 1 j a s . The dependence of Bi on a s is a func- 
tion of the dimensionless parameter ~ y/MBia s <C 1, 
and therefore take to be zero. 

In order to account for the theoretical uncertainty in 
the EFT calculation we assign an additional 10% random 
variation to the bindings of 3 He and 3 H (computed at 
NLO in EFT) and a 30% variation on the value of the 
4 He binding (computed at LO only), as will be shown 
more explicitly below. 

The reaction n + p <-> d + 7 was extensively analyzed 
in Ref. [37] using a N 4 LO calculation in the pionless 
EFT. The inputs at this order are the scattering length 
a s , the deuteron binding energy, the corresponding ef- 
fective range parameters, the magnetic moments of the 
deuteron, and a single two-nucleon-one-photon term fixed 
by experiment. We use the variation of B 2 and a s given 
in Eqs. ( 1.1 1 and (1 1 -2b to compute, with the help of the 



explicit formula in|27 , the relative change in the reac- 
tivity as a function of the temperature and use this as 
input for the BBN code. In [?] it was argued that the 
reactivity (av) scales as ~ B% a 2 . We verified with the 



dww\/w 2 — l(wo—w)' 



2na 



1 



Vw 2 - 1 1 _ g V /^TTT 

(1.9) 

which describes the phase space and the Coulomb repul- 
sion. The variation of T with the quark masses is given 
then by 



dT 



r dm„ 



—f(- 

dm V 



A 



m g 3 d(g A ) 
1 + 3g A dm q 



(1.10) 

The dependence of g A with the quark mass is given, at 
NLO in chiral perturbation theory, by [29] 



9a 



9 a 



1 



9g 2 A ml m n 



32rr 2 F 2 bg ^A 



(9 A ~ ^) m l 
32tt 2 F 2 



log(^f) 



m 

A 
(1.11) 

where g A is the chiral value of g A , F 93 MeV and A, A' 
are constants of order 1 GeV dependent on the Gasser- 
Leutwyler coefficients [3D]. Numerically we find 



1 + 3g A dm, 



1 

2T 



2 

A) 

3g\ drriji 



0.2. 



(1.12) 



The variation of the phase space f(A/m e ) with the 
quark masses can be estimated as 



m q /(^) 
/(£) dm i 



2/(£) dm v 
m-Tr df(w ) 
2f(w ) dw L °= 



(1.13) 



dA/m e 
dm^ 



The value of f(wo) and its derivative at wq = A/m e is 
found numerically to be 1.64 and 4.25, respectively. The 
variation of A/ m e with m q can be estimated by splitting 
A into a strong interaction component A s proportional 
to the up and down quark mass difference (and, con- 
sequently, to the value of m q ) and an electromagnetic 
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piece A e . m ., largely independent of m q . Unfortunately, 
the electromagnetic part is due to short distance effects 
and cannot be directly computed in a reliable way. The 
best handle we have on its value comes from chiral pertur- 
bation theory, where the up and down quark mass ratio, 
the meson spectrum, and the best estimate of the nucleon 
a— term are used as inputs to extract A s . The value ob- 
tained for A s in this manner is consistent with that ob- 
tained from lattice QCD calculation pJTJ. The difference 
between this value of A s and the measured value of the 
neutron-proton mass splitting gives A e m — — 0.76±0.30 

rni- 



ii. RESULTS 



In order to deal with the highly non-linear dependence 
of the final abundances on the quark masses and, at same 
time, to include estimates of theoretical errors, we use a 
stochastic procedure. More specifically, for a given quark 
mass variation Sm q /m q , we specify the binding energies 
of 2 H , 3 H , 3 He , 4 He , the reactivity (cry) for n+p <-» d + 
7, the neutron lifetime r and the phase space parameter 
Q. All other BBN parameters are kept at their present 
values. 



Chiral perturbation theory predicts a quark mass de- 
pendence of A s of the form A s = Am„(md — m u )/(m<j + 
m u ), a formula valid up to NLO since the leading order 
loop contribution to the nucleon mass cancels between 
the neutron and proton. We then have 



1 df(w ) 
~ f(w ) dw lw °= 
_ 1 df(w ) 
~ f(w ) dw U= 



2m fi 
A, 



A 



m d - to. 



TO,, 



2to w 



10.4 ±1.5. (1.14) 



Notice that we are taking both the up and down mass 
to vary while keeping the ratio md/m u fixed. As the 
dependence in Eq. (1.14) dominates over Eq. 
finally find 



(1.121 



r dm„ 



10.6 ± 1.5. 



(1.15) 



The quark mass variation of the neutron lifetime is rel- 
evant for our calculation. In order to sec that, let us 
remember that the neutron number, after the weak in- 
teractions are decoupled, decreases until BBN starts at 
t « 168 s. The suppression factor in standard BBN 
is thus e - 168 / 885 w 0.827. A 5% in creas e of quark 
masses would lead, according to Eq. (1.14), to a de- 
crease, of about 50% in the neutron lifetime and the 
suppression factor would change to e~ 252 / 885 rj 0.752, 
leading to a 4 He abundance change of about 10% in the 
4 He abundance, a variation comparable to the observa- 
tional uncertainties. 



The rate of other weak reactions changes in a similar 
manner. The phase space integrals are more involved and 
are, in BBN codes, computed "on the fly", taking the 
ratio Q = A/m e as input. We calculated the variation 
of Q as 



ni q dQ 
Q dm q 



dA 



2A gJtott 
1.59 ±0.23 



A. 
A 



(1.16) 



We have randomly generated a set of 300 values of 
scattering lengths a s , deuteron bindings B 2 with a Gaus- 
sian distribution with mean value and standard deviation 
given by 



X = 



ox 



1 



m q dX 


m q 


dX 


\ Am 


X dm q 
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dm q 


) m q 


m q dX 


m q 


dX 


\ Am 


X dm q 


X 


dm q 


/ m q 



A phys , 

j^-phys ^ 

(2.1) 



where X stands for either a s or B 2 and the "+" and "-" 
subscripts refer to the higher and lower values of dX/dm q 
allowed by Eqs. ( 1.1 ) and ( 1.2 ). The variations of a s and 
B 2 are assumed to be uncorrelated. From the ensemble 
of a s and B 2 obtained as above, we compute a corre- 



sponding ensemble of binding energies using Eq. (1.5) 
and add to the result a 10% (for 3 H and 3 He ) or 30% 
(for 4 He ) relative random error in order to take into ac- 
count theoretical errors discussed in the previous section. 
The binding energies of 3 H , 3 He and 4 He are then given 
by 



a s dBi 
Bs, da. 



Bo dB, 



(a.s 



7 phys 



Bi dBo 

(2.2) 

where i indexes the three nuclei 3 H , 3 He and 4 He , the 
superscript "phys" stands for the present, experimental 
values of the quantity, are Gaussian random variables 
with central value and standard deviation equal to 1, 
and ti is the theoretical error of the extrapolation equal 
to 0.1 (for 3 H and 3 He ) and 0.3 (for 4 He ). 

Similarly, the reactivity (<tv)t of the n + p — » d + 7 
reaction was computed as a function of the temperature 
T for the ensemble of a s , B 2 values determined by Eq. |2.1| 
using the explicit expression for the cross section from 
[27] . The high order expansion of this calculation in [27] 
is accompanied with very small theoretical errors that we 
subsequently neglect. 

We also generated, for each value of 5m q /m q , a set of 
300 random values of r and Q whose distribution reflect 
the discussion in the previous section. More specifically, 
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FIG. 2: The yellow bands show the (1-cr allowed abundances 
for 4 He and 2 H . The two clouds show the result of 300 sim- 
ulations, both with Am q /m q = -1 % but two different values 
of 7)io. The lower cloud (ochre online) is the result of taking 
r/io — 6.23 and the upper cloud (burgundy online) the value 
7710 = 4.60. There is very little change on the 4 He yield but 
the deuterium yield changes enough to render the deuterium 
abundance useless in putting a constraint on Am q /m q . 



these values were generated through the formula 



1 



T phys 

np h y s 



1 + (10.6 + 1.5$)^ 
m q 

1 + (1.59 + 0.23£)^2 
m„ 



(2.3) 



where £ is a Gaussian random variable with central value 
and standard deviation 1. Notice that this £ is indepen- 
dent of the £j used in the determination of the binding 
energies but the same £ is used in both r and Q since the 
leading theoretical uncertainties on both quantities stem 
from the same determination of the cr-term. 

For a given value of 8m q /m q , a set of values for 
B2, B3 H , B3 He , B4, (ctv)t was paired one of set of r and Q 
values and used in a standard BBN code. The BBN code 
we have used in our analysis is based on Refs. [331 134] 
and is publicly available [35] • The code was modified to 
accept temperature-dependent variations in the reactiv- 
ity corresponding to the n + p — > d + 7 reaction and the 
rate of weak interaction processes was changed according 
to Eq. ( [L15{ and Eq. ( flTo) . The Q-values of all BBN 
reactions with 2 H , 3 H , d He , and 4 He as either parent 
or daughter products of reactions were allowed to vary in 
accordance with the changes in binding energies of these 
nuclei. The baryon-to-photon ratio 77 was changed over 
a range discussed below. Otherwise, the standard input 
parameters were used in our BBN simulations. 

The main feature seen in the simulations is that a vari- 
ation in 77 shifts the deuterium abundance but has little 
effect on the 4 He yields (see Fig. (|2|). A larger value of 77 
implies in a larger baryon density, a more complete burn- 
ing of the neutrons into 4 He nuclei and a smaller deu- 
terium abundance. As a consequence, in the absence of 
a restriction on the value of 77 from other considerations, 
the deuterium abundance does not put any constraints 



FIG. 3: The yellow bands show the (1-cr allowed abundances 
for 4 He and 2 H . The five clouds show the result of 300 
simulations at each one of the values (from left to right): 
Am q /m q = 2% (green online), Am q /m q — 0.7% (blue on- 
line), Am q /m q = —0.5% (purple online), Am q /m q = — 1% 
(ochre online), and Am q /m q — —2% (red online). 



on the range of allowed quark masses variations. 

Additional constraints on the value of 77 come from 
studies of the large-scale structure of the Universe. The 
actual numerical value of the constraints, however, de- 
pends on assumptions made in these analyses, includ- 
ing assumptions on the initial spectrum of fluctuations. 
For instance, the lower range of the determination of 
7710 = 4.79 ± 0.019 in [36] and the central value of the 
determination of r)io — 6.23 ± 0.17 in [37], are shown 
for 5m q /m q in Fig. A similar plot results for other 
values of 5m q /m q . Consequently, any reasonable change 
in the deuterium abundance can be accommodated by 
a change in the value of 7710. If we restrict ourselves to 
the much narrower range 7710 = 6.23 ± 0.17 [37j . the deu- 
terium abundances can play a role. However, the values 
in the range 7710 = 6.23 ± 0.17 are in tension with the 
observed deuterium abundances. BBN, by itself, prefers 
the slightly lower range 5.1 < 7710 < 6.5, at the 95% 
confidence level [38j. Thus, even with the current physi- 
cal values of m q , the predicted deuterium abundance lies 
just outside the 1 — a band, making it difficult to distin- 
guish the allowed and forbidden values of m q based on 
Yi. Therefore, to proceed further, we disregard the deu- 
terium abundances and look at how the 4 He abundances 
change with the quark masses. 

In Fig. ([3]) we show the result of changing the quark 
masses by five values: 2%, 0.7%, -0.5%, -1% and -2%, 
all corresponding to 7710 = 6.23. Each one of these val- 
ues of 5m q /m q is represented by a cloud of points in 
the Y4 x Yi plane. The spread between the 300 points 
in each cloud accounts for the theoretical uncertainties in 
the extrapolation of the parameter inputs as described by 
Eqs. (2.2 1 and (2.3). The tendency is for a smaller Y4 for 
Two main mechanisms account for 



larger values of m q . 
this general trend. First, large values of m q imply larger 
values of A, as well as a larger phase space for neutron de- 
cay and therefore shorter neutron lifetime. Consequently, 
more neutrons decay by the time BBN starts the assem- 
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bly of 4 He , resulting in smaller 4 He yields. In addition, 
Eq. (1.1 1 shows that a larger m q implies a smaller B^. 



The deuteron, being less bound, takes longer to form, de- 
laying the onset of 4 He formation and giving even more 
time for the neutrons to decay, reducing further the He 
yields. There is also a weak tendency to have smaller Y2 
for smaller m q , a trend not so easily explained. 

Based on the data shown on Fig. ^ we put a bound 
on the allowed values of quark mass changes at 



1% < ^ < 0.7%, 



(2.4) 



which is the main result of this paper. We refrain from 
assigning a numerical value to the uncertainty in this es- 
timate, as an attempt in this direction would require us 
to assign a precise statistical meaning to our theoretical 
uncertainties. While there are reasons to take these un- 
certainties seriously at the qualitative level, we believe 
them to be superior to the model calculations used pre- 
viously. 



III. CONCLUSION 

We have estimated the abundances of 2 H and 4 He pro- 
duced in the standard BBN scenario under the assump- 
tion that the light quark masses were shifted at the BBN 
time from their present values. In order to perform this 
calculation we have used input from several effective field 
theories as well as lattice QCD results to connect the 
quark mass variation to the relevant nuclear physics per- 
tinent to BBN. We found that a variation beyond the 
— 1% ^ ~ 0-7% range to be likely inconsistent with 
the observed abundances. 

Two of the BBN parameters played the largest roles in 
changing the light element yields: the deuteron binding 
energy B2 (with the 3 H , 3 He and 4 He binding energies 
strongly correlated with B2) and the neutron lifetime. 
The dependence of the neutron lifetime on the quark 
mass values is well constrained by theory. The varia- 
tion of the deuteron binding is, however, much less con- 
strained and several venues of further progress are clearly 
visible (for a very recent study, see 39J ) . Lattice QCD 
calculations of nucleon-nucleon interactions, even if per- 
formed at unphysical values of m q , would go a long way 
in narrowing these constraints. As long as they are per- 
formed with quark masses low enough to be within the 
region of validity of the chiral nuclear EFT, they deter- 
mine reliably the value of parameters of the EFT neces- 
sary for the extrapolation of the deuteron binding energy. 
The binding energies of 3 H , 3 He and, specially 4 He , 
can and should be computed in the pionless nuclear EFT 
to higher orders so that theoretical uncertainties associ- 



ated with these quantities decrease. Finally, a better un- 
derstanding of the quark mass variation of other thresh- 
old parameters like effective ranges, magnetic moments, 
etc. . . , would also allow for a more precisely constrained 
calculation of the binding energies on nuclei larger then 
4 He. 

A number of other works have also considered the ef- 
fects of a variation of the quark masses on properties 
of light nuclei. For example, in Ref. [5] this effect was 
implemented by a change of the pion mass in the phe- 
nomenological model interaction employed in the calcu- 
lations. Such an interaction is only able to capture the 
true quark mass dependence to a limited degree since it 
is not constructed as a systematic expansion in powers of 
I Kqqd- In particular, quark mass dependent short- 
distance contact operators (such as the D-term) discussed 
in the above text do not appear in standard phenomeno- 
logical interactions. 

Since we are not presently able to obtain reliable val- 
ues for the 7 Li binding energies, the 7 Li abundances we 
compute are not very meaningful and were not used to 
put constraints on the quark mass variations. Future ad- 
vances in the nuclear pionless effective theory may change 
this and allow us to address the "Lithium problem" as a 
signal of quark mass variation. 

Finally, it should be pointed out that in models of 
physics beyond the standard model, the value of the 
quark masses are derived quantities and variations of 
them may well be correlated with other quantities. In 
particular, it may seem unnatural to expect the masses 
of different quark flavors to vary together, unless this 
variation is being driven by a change in the Higgs vac- 
uum expectation value. If that is the change in 
the quark masses will be correlated with changes in the 
vector boson masses, changing the strength of strong in- 
teraction at low energies. The effect of those changes on 
BBN can be easily tracked in a manner similar to what 
was done in this paper. We plan to consider BBN bounds 
on the Higgs vacuum expectation value change in a future 
publication. 
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